#include <iostream>
#include <sstream>
#include <math.h>
#include <vector>
#include <string>

#include <TVirtualMC.h>
#include <TGeoManager.h>
#include <TGeoShape.h>
#include <TGeoBBox.h>
#include <TGeoTube.h>
#include <TGeoMatrix.h>

#include "sPHENIX_DetectorConstruction.h"

#define PI 3.14159265

/// \cond CLASSIMP
ClassImp(sPHENIX_DetectorConstruction)
/// \endcond

using namespace std;

//_____________________________________________________________________________
sPHENIX_DetectorConstruction::sPHENIX_DetectorConstruction()
  : TObject(),
    fWorldWidth(0.), fWorldLength(0.),  
    fEMCalAbsorberLength(0.), fEMCalSensingLength(0.), fEMCalTowerWidth(0.),
    EMCalNofCells(0.), EMCalNofTowersAxis(0.), EMCalNofTowersHole(0.), EMCalLocation(0.),
    fBeamROuter(0.),
    NofGEMs(0.),fGEMLength(0.), fGEMThickness(0.), fGEMRapidity(0.),
    GEMLocation(0.),
    PreshowerLocation(0.),
    fHCalAbsorberLength(0.), fHCalSensingLength(0.), fHCalTowerWidth(0.),
    HCalNofCells(0.), HCalNofTowersAxis(0.), HCalNofTowersHole(0.), HCalLocation(0.),
  ImedAir(0),ImedPb(0),ImedSc(0),ImedC(0),ImedAr(0), ImedFe(0)
{
  /// Default constuctor

  //fpMagField = new ExN02MagneticField();

  //--------- Sizes of the principal geometrical components (solids)  ---------
 
  // - World  
  fWorldWidth = 800;
  fWorldLength = 1200.;

  // - EM Calorimeter
  fEMCalAbsorberLength = 0.15;
  fEMCalSensingLength = 0.4; 
  fEMCalTowerWidth = 5.535;      // [cm]
  EMCalNofCells = 66;
  EMCalNofTowersAxis = 126;      // number
  EMCalNofTowersHole = 16;
  EMCalLocation = 300;   //300

  // - Beam Pipe
  fBeamROuter = 3.;

  // - GEMS
  NofGEMs = 4;
  fGEMLength = 1.2;   
  fGEMThickness = 0.1;
  fGEMRapidity = 1.;             // rapidity
  GEMLocation.push_back(70.);   // GEM0 location
  GEMLocation.push_back(80.);   // GEM1 location
  GEMLocation.push_back(140.);  // GEM2 location
  GEMLocation.push_back(150.);  // GEM3 location

  // - Hadronic Calorimeter
  fHCalAbsorberLength = 2.5;
  fHCalSensingLength = 0.5; 
  fHCalTowerWidth = 11.07;      // [cm]
  HCalNofCells = 20;
  HCalNofTowersAxis = 64;      // number
  HCalNofTowersHole = 8;
  HCalLocation = 360;       // 360

  // - Preshower
  PreshowerLocation = 297.;
  fPreshowerLength = 1.;
}

//_____________________________________________________________________________
sPHENIX_DetectorConstruction::~sPHENIX_DetectorConstruction()
{
  /// Destructor

  // delete fpMagField;
}

//_____________________________________________________________________________
void sPHENIX_DetectorConstruction::ConstructMaterials()
{
  /// Construct materials using TGeo modeller

  //--------- Material definition ---------

  // Create Root geometry manager 
  new TGeoManager("sPHENIX_geometry", "E02 VMC example geometry");
   
  Double_t a;        // Mass of a mole in g/mole   
  Double_t z;        // Atomic number
  Double_t density;  // Material density in g/cm3
 
  // Elements
  
  TGeoElement* elN  = new TGeoElement("Nitrogen", "N", z= 7., a= 14.01);
  TGeoElement* elO  = new TGeoElement("Oxygen"  , "O", z= 8., a= 16.00);

  // Materials

  TGeoMixture* matAir
    = new TGeoMixture("Air", 2, density = 1.29e-03);
  matAir->AddElement(elN, 0.7); 
  matAir->AddElement(elO, 0.3); 

  TGeoMaterial* matLead 
    = new TGeoMaterial("Lead", a = 207.19, z = 82., density = 11.35); 

  TGeoMaterial* matSc             
    = new TGeoMaterial("SensingMat", a = 6.67, z = 3.6, density = 1.18);
 
  TGeoMaterial* matBe
    = new TGeoMaterial("Beryllium", a = 9.01, z = 4., density = 1.85);   

  TGeoMaterial* matC
    = new TGeoMaterial("Carbon", a = 12.01, z = 6., density = 2.267);   

  TGeoMaterial* matAr
    = new TGeoMaterial("Argon", a = 39.95, z = 18., density = 1.67);   

  TGeoMaterial* matFe
    = new TGeoMaterial("Iron", a = 55.85, z = 26., density = 7.87);   

  // Tracking media

  Double_t param[20];
  param[0] = 0;     // isvol  - Not used
  param[1] = 2;     // ifield - User defined magnetic field
  param[2] = 10.;   // fieldm - Maximum field value (in kiloGauss)
  param[3] = -20.;  // tmaxfd - Maximum angle due to field deflection 
  param[4] = -0.01; // stemax - Maximum displacement for multiple scat 
  param[5] = -.3;   // deemax - Maximum fractional energy loss, DLS 
  param[6] = .001;  // epsil - Tracking precision
  param[7] = -.8;   // stmin
  for ( Int_t i=8; i<20; ++i) param[i] = 0.;

  ImedAir = 1; 
  new TGeoMedium("Air", ImedAir, matAir, param);
  
  ImedPb = 2;
  new TGeoMedium("Lead", ImedPb, matLead, param);
  
  ImedSc = 3;
  new TGeoMedium("SensingMat", ImedSc, matSc, param);

  ImedBe = 4;
  new TGeoMedium("Beryllium", ImedBe, matBe, param);

  ImedC = 5;
  new TGeoMedium("Carbon", ImedC, matC, param);
  
  ImedAr = 6;
  new TGeoMedium("Argon", ImedAr, matAr, param);

  ImedFe = 7;
  new TGeoMedium("Iron", ImedFe, matFe, param);
}    

//_____________________________________________________________________________
void sPHENIX_DetectorConstruction::ConstructGeometry()
{
  /// Contruct volumes using TGeo modeller

  //------------------------------ 
  // World
  //------------------------------
 
  Double_t world[3] = {fWorldWidth/2,fWorldWidth/2,fWorldLength/2};
  TGeoVolume *top = gGeoManager->Volume("WORLD","BOX",ImedAir,world,3);
  gGeoManager->SetTopVolume(top);

  //TGeoVolume *magnetvol = new TGeoVolume("magvol",magnet);
  //top->AddNode(magnetvol,1,new TGeoTranslation(0,0,0));

  //------------------------------ 
  // Cell - Node(s) of EM Tower
  //------------------------------ 

  Double_t fEMCalCellLength = fEMCalSensingLength+fEMCalAbsorberLength;  ///< Full length of cell
  
  TGeoBBox *emcalcell = new TGeoBBox( 0.5*fEMCalTowerWidth, 
				      0.5*fEMCalTowerWidth,
				      0.5*fEMCalCellLength);
  TGeoVolume *emcalcellvol = new TGeoVolume("emcalcellvol",emcalcell);
  emcalcellvol->SetMedium(gGeoManager->GetMedium(ImedAir));


  //------------------------------ 
  // Absorber - Node of EM Cell
  //------------------------------ 

  TGeoBBox *emcalabsorber = new TGeoBBox( 0.5*fEMCalTowerWidth, 
					  0.5*fEMCalTowerWidth,
					  0.5*fEMCalAbsorberLength);
  TGeoVolume *emcalabsorbervol = new TGeoVolume("emcalabsorbervol",emcalabsorber);
  emcalabsorbervol->SetMedium(gGeoManager->GetMedium(ImedPb));
  emcalcellvol->AddNode(emcalabsorbervol,1,new TGeoTranslation(0.,0.,-0.5*fEMCalCellLength + 0.5*fEMCalAbsorberLength));

  //------------------------------ 
  // SensingMat - Node of EM Cell
  //------------------------------ 

  TGeoBBox *emcalsens = new TGeoBBox( 0.5*fEMCalTowerWidth, 
				      0.5*fEMCalTowerWidth,
				      0.5*fEMCalSensingLength);
  TGeoVolume *emcalsensvol = new TGeoVolume("emcalsensvol",emcalsens);
  emcalsensvol->SetMedium(gGeoManager->GetMedium(ImedSc));
  emcalcellvol->AddNode(emcalsensvol,1,
			new TGeoTranslation(0,0,0.5*fEMCalCellLength - 0.5*fEMCalSensingLength));

  //------------------------------ 
  // Tower - Node(s) of EM Calorimeter
  //------------------------------ 
  
  Double_t fEMCalTowerLength = EMCalNofCells*fEMCalCellLength;      ///< Full length of tower

  TGeoBBox *emcaltower = new TGeoBBox( 0.5*fEMCalTowerWidth, 
				       0.5*fEMCalTowerWidth,
				       0.5*fEMCalTowerLength);
  TGeoVolume *emcaltowervol = new TGeoVolume("emcaltowervol",emcaltower);
  emcaltowervol->SetMedium(gGeoManager->GetMedium(ImedAir));  

  // Cell Placement in tower

  for( Int_t i = 0; i < EMCalNofCells; i++)
    emcaltowervol->AddNode(emcalcellvol,i,
			   new TGeoTranslation(0,0,-0.5*fEMCalTowerLength + (0.5 + i)*fEMCalCellLength));


  //------------------------------ 
  // EM Calorimeter
  //------------------------------ 

  TGeoBBox *emcal = new TGeoBBox( 0.5*fEMCalTowerWidth*EMCalNofTowersAxis, 
				  0.5*fEMCalTowerWidth*EMCalNofTowersAxis,
				  0.5*fEMCalTowerLength);
  TGeoVolume *emcalvol = new TGeoVolume("emcalvol",emcal);
  emcalvol->SetMedium(gGeoManager->GetMedium(ImedAir));


  // Tower Placement in calorimeter

  Double_t coordstart = -(EMCalNofTowersAxis/2*fEMCalTowerWidth) + 0.5*fEMCalTowerWidth;
  Double_t xc, yc;
  
  Int_t EMCalNofTowerHoleLower = EMCalNofTowersAxis/2 - EMCalNofTowersHole/2;
  Int_t EMCalNofTowerHoleUpper = EMCalNofTowerHoleLower + EMCalNofTowersHole -1; 
  Int_t counter = 0;

  for(Int_t i = 0; i < EMCalNofTowersAxis; i++){
    yc = coordstart + i*fEMCalTowerWidth;
    for(Int_t j = 0; j < EMCalNofTowersAxis; j++){
      xc = coordstart + j*fEMCalTowerWidth;
      if( !((i >= EMCalNofTowerHoleLower && i <= EMCalNofTowerHoleUpper) 
	    && (j >= EMCalNofTowerHoleLower && j <= EMCalNofTowerHoleUpper)))
	emcalvol->AddNode(emcaltowervol,counter,new TGeoTranslation(xc,yc,0.));
      counter++;
    }
  }

  top->AddNode(emcalvol,1, new TGeoTranslation(0., 0., EMCalLocation + fEMCalTowerLength/2));

  emcaltowervol->SetVisibility(kTRUE);

  //------------------------------ 
  // GEMS
  //------------------------------ 

  Double_t fGEMRInner = fBeamROuter;      // For confusion reasons
  Double_t fGEMAngle = 
    atan(exp(-fGEMRapidity))*2;           // Radians
  Double_t fGEMGasLength = 
    fGEMLength - 2*fGEMThickness;         // Full GEMGas Length
  Double_t plateZ = -fGEMLength/2+fGEMThickness/2;    // Location of Plate inside mother

  vector<Double_t>fGEMROuter (NofGEMs);
  vector<TGeoTube*> gems (NofGEMs);
  vector<TGeoVolume*> gemvols (NofGEMs);
  vector<TGeoTube*> fsgems (NofGEMs);
  vector<TGeoVolume*> fsgemvols (NofGEMs);
  vector<TGeoTube*> gemgass (NofGEMs);
  vector<TGeoVolume*> gemgasvols (NofGEMs);

  for(int i = 0; i < NofGEMs; i++){
    
    // --- GEM
    fGEMROuter[i] = GEMLocation[i] * tan(fGEMAngle);
    gems[i] = new TGeoTube( fGEMRInner,   // RInner
			    fGEMROuter[i],
			    0.5*fGEMLength);
    gemvols[i] = new TGeoVolume("gemvol",gems[i]);
    gemvols[i]->SetMedium(gGeoManager->GetMedium(ImedC));
    
    // --- Sensative Front
    fGEMROuter[i] = GEMLocation[i] * tan(fGEMAngle);
    fsgems[i] = new TGeoTube( fGEMRInner + fGEMThickness,   // RInner
			      fGEMROuter[i] - fGEMThickness,
			      0.5*fGEMThickness);
    fsgemvols[i] = new TGeoVolume("fsgemvol",gems[i]);
    fsgemvols[i]->SetMedium(gGeoManager->GetMedium(ImedC));
    
    gemvols[i]->AddNode(fsgemvols[i],i,new TGeoTranslation(0., 0., plateZ));


    // --- Gas
    gemgass[i] = new TGeoTube( fGEMRInner + fGEMThickness,   // RInner
			       fGEMROuter[i] - fGEMThickness,
			       0.5*fGEMGasLength);
    gemgasvols[i] = new TGeoVolume("gemgasvol",gemgass[i]);
    gemgasvols[i]->SetMedium(gGeoManager->GetMedium(ImedAr));
    
    gemvols[i]->AddNode(gemgasvols[i],i,new TGeoTranslation(0., 0., 0.));
    
    top->AddNode(gemvols[i],i, new TGeoTranslation(0., 0., 
						   GEMLocation[i]+fGEMLength/2));
  }

  //------------------------------ 
  // Preshower Tower
  //------------------------------  

  Double_t PreshowerWidth = fEMCalTowerWidth;     // Not to get confused
  
  TGeoBBox *preshowertower = new TGeoBBox( 0.5*PreshowerWidth, 
					   0.5*PreshowerWidth,
					   0.5*fPreshowerLength);
  TGeoVolume *preshowertowervol = new TGeoVolume("preshowertowervol",preshowertower);
  preshowertowervol->SetMedium(gGeoManager->GetMedium(ImedSc));


  //------------------------------ 
  // Preshower
  //------------------------------ 

  Int_t PreshowerNofTowersAxis = EMCalNofTowersAxis;

  TGeoBBox *preshower = new TGeoBBox( 0.5*PreshowerWidth*PreshowerNofTowersAxis, 
				      0.5*PreshowerWidth*PreshowerNofTowersAxis,
				      0.5*fPreshowerLength);
  TGeoVolume *preshowervol = new TGeoVolume("preshowervol",preshower);
  preshowervol->SetMedium(gGeoManager->GetMedium(ImedAir));

  // Preshower Placement

  coordstart = -(PreshowerNofTowersAxis/2*PreshowerWidth) + 0.5*PreshowerWidth;

  Int_t fNofPreshowersHole = EMCalNofTowersHole;
  Int_t fNofPreshowerHoleLower = PreshowerNofTowersAxis/2 - fNofPreshowersHole/2;
  Int_t fNofPreshowerHoleUpper = fNofPreshowerHoleLower + fNofPreshowersHole -1; 
  counter = 0;

  for(Int_t i = 0; i < PreshowerNofTowersAxis; i++){
    yc = coordstart + i*PreshowerWidth;
    for(Int_t j = 0; j < PreshowerNofTowersAxis; j++){
      xc = coordstart + j*PreshowerWidth;
      if( !((i >= fNofPreshowerHoleLower && i <= fNofPreshowerHoleUpper) 
	    && (j >= fNofPreshowerHoleLower && j <= fNofPreshowerHoleUpper)))
	preshowervol->AddNode(preshowertowervol,counter,new TGeoTranslation(xc,yc,0.));
      counter++;
    }
  }

  top->AddNode(preshowervol,1, new TGeoTranslation(0., 0., PreshowerLocation + fPreshowerLength/2));

  //------------------------------ 
  // Cell - Node(s) of Hardonic Tower
  //------------------------------ 

  Double_t fHCalCellLength = fHCalSensingLength+fHCalAbsorberLength;  ///< Full length of cell
  
  TGeoBBox *hcalcell = new TGeoBBox( 0.5*fHCalTowerWidth, 
				      0.5*fHCalTowerWidth,
				      0.5*fHCalCellLength);
  TGeoVolume *hcalcellvol = new TGeoVolume("hcalcellvol",hcalcell);
  hcalcellvol->SetMedium(gGeoManager->GetMedium(ImedAir));


  //------------------------------ 
  // Absorber - Node of Hadronic Cell
  //------------------------------ 

  TGeoBBox *hcalabsorber = new TGeoBBox( 0.5*fHCalTowerWidth, 
					  0.5*fHCalTowerWidth,
					  0.5*fHCalAbsorberLength);
  TGeoVolume *hcalabsorbervol = new TGeoVolume("hcalabsorbervol",hcalabsorber);
  hcalabsorbervol->SetMedium(gGeoManager->GetMedium(ImedFe));
  hcalcellvol->AddNode(hcalabsorbervol,1,new TGeoTranslation(0.,0.,-0.5*fHCalCellLength + 0.5*fHCalAbsorberLength));

  //------------------------------ 
  // SensingMat - Node of Hadronic Cell
  //------------------------------ 

  TGeoBBox *hcalsens = new TGeoBBox( 0.5*fHCalTowerWidth, 
				      0.5*fHCalTowerWidth,
				      0.5*fHCalSensingLength);
  TGeoVolume *hcalsensvol = new TGeoVolume("hcalsensvol",hcalsens);
  hcalsensvol->SetMedium(gGeoManager->GetMedium(ImedSc));
  hcalcellvol->AddNode(hcalsensvol,1,
			new TGeoTranslation(0,0,0.5*fHCalCellLength - 0.5*fHCalSensingLength));

  //------------------------------ 
  // Tower - Node(s) of Hadronic Calorimeter
  //------------------------------ 
  
  Double_t fHCalTowerLength = HCalNofCells*fHCalCellLength;      ///< Full length of tower

  TGeoBBox *hcaltower = new TGeoBBox( 0.5*fHCalTowerWidth, 
				       0.5*fHCalTowerWidth,
				       0.5*fHCalTowerLength);
  TGeoVolume *hcaltowervol = new TGeoVolume("hcaltowervol",hcaltower);
  hcaltowervol->SetMedium(gGeoManager->GetMedium(ImedAir));  

  // Cell Plachent in tower

  for( Int_t i = 0; i < HCalNofCells; i++)
    hcaltowervol->AddNode(hcalcellvol,i,
			   new TGeoTranslation(0,0,-0.5*fHCalTowerLength + (0.5 + i)*fHCalCellLength));


  //------------------------------ 
  // Hadronic Calorimeter
  //------------------------------ 


  TGeoBBox *hcal = new TGeoBBox( 0.5*fHCalTowerWidth*HCalNofTowersAxis, 
				  0.5*fHCalTowerWidth*HCalNofTowersAxis,
				  0.5*fHCalTowerLength);
  TGeoVolume *hcalvol = new TGeoVolume("hcalvol",hcal);
  hcalvol->SetMedium(gGeoManager->GetMedium(ImedAir));

  // Tower Plachent in calorimeter

  coordstart = -(HCalNofTowersAxis/2*fHCalTowerWidth) + 0.5*fHCalTowerWidth;
  xc, yc;
  
  Int_t HCalNofTowerHoleLower = HCalNofTowersAxis/2 - HCalNofTowersHole/2;
  Int_t HCalNofTowerHoleUpper = HCalNofTowerHoleLower + HCalNofTowersHole -1; 
  counter = 0;

  for(Int_t i = 0; i < HCalNofTowersAxis; i++){
    yc = coordstart + i*fHCalTowerWidth;
    for(Int_t j = 0; j < HCalNofTowersAxis; j++){
      xc = coordstart + j*fHCalTowerWidth;
      if( !((i >= HCalNofTowerHoleLower && i <= HCalNofTowerHoleUpper) 
	    && (j >= HCalNofTowerHoleLower && j <= HCalNofTowerHoleUpper)))
	hcalvol->AddNode(hcaltowervol,counter,new TGeoTranslation(xc,yc,0.));
      counter++;
    }
  }

  top->AddNode(hcalvol,1, new TGeoTranslation(0., 0., HCalLocation + fHCalTowerLength/2));

  //------------------------------ 
  // Beam Pipe
  //------------------------------ 

  TGeoTube *beampipe = new TGeoTube( 2.0,   // RInner
				     fBeamROuter,
				     0.5*fWorldLength);
  TGeoVolume *beampipevol = new TGeoVolume("beampipevol",beampipe);
  beampipevol->SetMedium(gGeoManager->GetMedium(ImedBe));

  top->AddNode(beampipevol,1,new TGeoTranslation(0., 0., 0.));


  // close geometry
  gGeoManager->CloseGeometry();
    
  // notify VMC about Root geometry
  gMC->SetRootGeometry();
}
//_____________________________________________________________________________
string sPHENIX_DetectorConstruction::ConvertNum(Int_t ii){
  stringstream ss;         //create a stringstream
  ss << ii;            //add number to the stream
  return ss.str();//return a string with the contents of the stream
}

